"""
    在执行提取操作前要先对与之无关的影像区域进行遮罩。这种技术被称为阈值分割（二值化）。本例中，
阈值图像上岛屿和背景区域的水域形成了足够鲜明的对比，阈值化（二值化）后能够完美地将它分隔开来。
"""
from osgeo import gdal, gdal_array, ogr, osr
import shapefile
import pngcanvas
# 地理坐标转换为屏幕坐标
def world2screen(bbox, w, h, x, y):
    minx, miny, maxx, maxy = bbox
    xdist = maxx - minx
    print(xdist,"max",maxx,"min",minx)
    ydist = maxy - miny
    xratio = w/xdist
    yratio = h/ydist
    px = int(w - ((maxx - x) * xratio))
    py = int((maxy - y) * yratio)  # 为了适应pngcanvas的坐标系，对y做了翻转
    return px, py

src = "TIF/islands/islands_classified.tiff"
output = "TIF/islands/extract.shp"
output_layer = "extract"
#
srcImage = gdal.Open(src)
band = srcImage.GetRasterBand(1)
mask = band
# geoTrans = srcImage.GetProjectionRef()
# print(geoTrans)
# 创建输出的Shp文件
driver = ogr.GetDriverByName("ESRI Shapefile")
shp = driver.CreateDataSource(output)
# 拷贝空间索引
srs = osr.SpatialReference()
srs.ImportFromWkt(srcImage.GetProjectionRef())
layer = shp.CreateLayer(output_layer, srs=srs)
# 创建dbf文件
fd = ogr.FieldDefn("DN", ogr.OFTInteger)
layer.CreateField(fd)
dst_field = 0
extract = gdal.Polygonize(band, mask, layer, dst_field, [], None)

width = 800
height = 600
r = shapefile.Reader(output)
print(r)
print(r.fields)
print(r.bbox)
polygons = []
for shape in r.shapes():
    for i in range(len(shape.parts)):
        pixels = []
        pt = None
        if i < len(shape.parts) - 1:
            pt = shape.points[shape.points[i]:shape.parts[i+1]]
        else:
            pt = shape.points[shape.parts[i]:]
        for x, y in pt:
            px, py = world2screen(r.bbox, width, height, x, y)
            pixels.append([px, py])
        polygons.append(pixels)
c = pngcanvas.PNGCanvas(width, height)
for p in polygons:
    c.polyline(p)
f = open("extract.png", "wb")
f.write(c.dump())
f.close()

r.close()